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10.0 OPERATION OF C0DUCT CODE 


10.1 Runstream 

The following runstream for a UNIVAC 1100 operating system is used to assign 
input/output disc files and to execute the C0DUCT coordinate generator code. 


@ASG,A 

FILE9. ,D/0/TRK/300000 

0USE 

9. ,FILE9. 

@ASG,A 

FILElO. .D/0/TRK/300000 

@USE 

10. .FILElO. 

@XQT 

C0DUCT 

« 


• 

data cards 

« 


@FIN 


FILE9 contains the coordinate data for a uniform mesh and FILElO contains 
the data for a nonuniform mesh. 
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10.2 Input Format 


The input format for the C0DUCT code is described on the input data coding 
forms which follow. With the exception of the first card (Title card) and the duct 
geometry cards, the input data cards follow in sets of three cards. The first of 
three is a blank separator card. The second of three is the input variable name 
and the third of three is the value of the input variable. In general the input 
data is read as follows : 

Card 1 Title Card 

Cards 2-4 Program Control Parameters 

Cards 5-7 Program Control Parameters 

Cards 8-10 Coordinate Generator Parameters 

Cards 11-13 Coordinate Generator Parameters 

Card 14 Number of duct geometry coordinates 

Cards 15 + Duct geometry coordinates 
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C0DUCT CODE INPUT 


Card 1 Title Card (18A4) 



I* 
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C0DUCT CODE INPUT 


Cards 2-4 Program Control Parameters (T6, 13, T20, II, T32, II, T43, 12, T53, F5.2, T64, E8.3) 



MAXIT Maximum number of iterations for conformal mapping solution 

IPRINT Print option: 

= 0 Print only input and final mapping 
= 1 Print iteration results 
= 2 Print approximate solution 

ICORD Coordinate generation option: 

= 0 Calculate conformal mapping only 
= 1 Calculate coordinates 

NIPOT ’ Approximate solution option 

= 0 Program determines number of approximate potential lines 
> 0 User specifies number of approximate potential lines 

DM Step control in approximate solution calculations 


ECONV 


Convergence criterion on conformal mapping iteration 
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C0DUCT CODE INPUT 


Cards 5-7 Program Control Parameters (T8, II, T18, 13, T30, 13) 



lESTIM Approximate solution option 

= 0 Calculate approximate solution only 
= 1 Calculate conformal mapping 

KN Number of streamlines to generate on uniform grid, number of integration steps 

NSD Potential line integration step control: 

=0 KN integration steps are used 

> 0 There will be NSD additional steps taken for each of the KN steps yielding 
NNS = (KN - 1) *NSD + KN total steps 
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C0DUCT CODE INPUT 


Cards 8-10 Coordinate Generation Control Parameters (T6, 13, T18, 13, T30, 13, T44, II, T51, E8.3, T68, II) 



JL Number of streamwise stations to output 

JLPTS Number of data points to output from smooth 
KL Number of nonuniformly spaced streamlines 

IGRID Grid option 

= 0 KN uniformly spaced streamlines output to unit lUUNIT 
=2 KL nonuniformly spaced streamlines output to unit INUNIT 
= 1 Both meshes (0) , (2) output 

LOP Control option for mesh distortion 

DDS Ratio of uniform/nonuniform mesh size at wall 


C(^DUCT CODE INPUT 


Cards 11-13 Coordinate Generation Control Parameters (T7, 12, T19, 12, T32, II, T43, 12, T55, II) 
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C0DUCT CODE INPUT 


Cards 15 -> Duct Geometry Coordinates (8F10.6) 




YU Upper wall y-coordinates 

XL Lower wall x-coordlnates 

YL Lower wall y-coordinates 

Where the data is read as follows 

READ (IRUNIT,40) (XU(J), J=1,NLF) 
READ (IRUNIT.40) (YU(J), J=1,NLF) 
READ (IRUNIT,40) (XL(J) , J=1,NLF) 
READ (IRUNIT,40) (Y.L(J) , J=1,NLF) 


40 FORMAT (8F10.6) 



10.3 Output Format 


The printed output from the C0DUCT code is given on the following pages and is 
largely self explanatory. These pages contain the names of the subroutines which 
calculate the data as well as any print options which may be involved. 
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Input Data Echo Page (1) 


Printed by RCNTRL 

Calculated by RCNTRL 

Options None 

Description 

Subroutine RCNTRL reads the input control data (cards 1-13) and prints the 
input control data. 
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Printed by 
Calculated by 
Options 


CORINP 

N/A 

None 


Description 

Subroutine CORINP reads the duct wall coordinates as described in Section 10.2 
and prints the wall coordinates. 


Pa^es 1-2 ^OTUT WALL_C00^^NATES 

Heading Variable Description 


PT# 

XU 

YU 

XL 

YL 


J 

^U 


Wall coordinate number 
Upper wall input X coordinate 
Upper wall input Y coordinate 
Lower wall input X coordinate 
Lower wall input Y coordinate 
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Pa£6£ 1-2 ^MOOTHED_W^L COORDINATES 

Heading Variable Description 


PT# 

XU 

YU 

XL 

YL 


Smoothed wall coordinate number 
Smoothed upper wall X coordinate 

Smoothed upper wall Y coordinate 

Smoothed lower wall X coordinate 

Smoothed lower wall Y coordinate 
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Subroutine MDAVIS determines the orientation of the duct and forces the inlet 
walls parallel. The rotational and scaling constant M and the exit wall angle are 
calculated and output by subroutine ROTATM. 


Paj^e_l £ALCULATED_DUCT GEOMETRY_PAR^IETEI« 

Heading Variable Description 

None Angle of inlet lower wall to horizon- 

tal 

None h True height of duct - perpendicular 

distance from lower wall to upper wall 

None M Rotational and scaling constant used 

in Schwartz-Cristoffel mapping 

None a„ Relative angle at exit of duct + r 

ri 
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Approximate Potential Flow Page (6) 


Printed by WALLV, ESTIMP 

Calculated by ESTIMP, ESTCOR, WALLV 

Options IPRINT 

IPRINT <2 no output 

IPRINT = 2 output 

Description 

Subroutivre ESTIMP solves for a geometric approximate potential solution in the 
duct using subroutine ESTCOR to determine the positions of a set of approximate po- 
tential lines. The wall velocities and curvatures at the end points of these lines 
are calculated and printed by subroutine WALLV. The solution along each wall is 
interpolated to the set of wall coordinates and an estimate of the Schwartz - 
Cristoffel mapping is determined. This estimate is calculated and printed by sub- 
routine ESTIMP. 


^a^e_l 


LOWER WALL ESTCOR FLOW ESTIMATE 


Heading 

Variable 

PT# 

NP 

ARCL 

S 

VMEAN 

V 

VWALL 

II 

KMEAN 

K 


Description 

Approximate potential line number 
Arc Length 

Duct midline velocity 
Lower wall velocity 
Duct midline curvature 
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Approximate Potential Flow Pages (Cont'd) 


Page 2 

UPPER WALL ESTCOR 

FLOW ESTIMATE 

Heading 

Variable 

Description 

PT# 

NP 

Approximate potential line number 

ARCL 

S 

Arc length 

VMEAN 

V 

Duct mid line velocity 

WALL 

Vu 

Upper wall velocity 

KMEAN 

K 

Duct mid line curvature 


Pages 3-A 

LOWER WALL 

POTENTIAL FLOW ESTIMATE 

Heading 

Variable 

Description 

PT# 

J 

Wall coordinate number 

WALL 

\ 

Lower wall velocity 

KWALL 

^L 

Lower wall curvature 

TT 

t 

t-plane pole estimates 

B 

b 

^-plane pole estimates 


Pa^es 5-6 

UPPER WALL 

POTENTIAL FLOW ESTIMATE 

Heading 

Variable 

Description 

PT# 

J 

Wall coordinate number 

WALL 

V 

Upper wall velocity 

KWALL 

Ku 

Upper wall curvature 

TT 

t 

t-plane pole estimate 

B 

b 

^-plane pole estimate 
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Iteration History Pages 


Printed by MDAVIS 

Calculated by MDAVIS, STEP 

Options IPRINT 

IPRINT = 0 not printed 

IPRINT > 0 history printed 

Description 

The current estimate of dz/dt is integrated along the duct walls fn subroutine 
MDAVIS by steps - tj in the t-plane. The results of the integration and the 

error of the current iteration are printed. 


Pa^e_l ITEMnON_H^STORY 


Heading 


Variable 


Description 


PT# 

TX 

TY 

X 

Y 

RATIO 

ERROR 

DZDT 


J 

t. 


R 


dz 

dt 


Wall coordinate number - j 
Real component of tj 
Imaginary component of tj 
Real component of Zj 
Imaginary component of Zj 
Arc length ratio 
Absolute error Zcj - Zj 
Complex derivative ot mapping 
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Iteration Summary Page (1) 


Printed by MDAVIS, CLOSUR 

Calculated by MDAVIS, CLOSUR 

Options None 

Description 

Subroutine MDAVIS calculates the maximum relative error in the coordinate cal- 
culation for each iteration and tests for convergence of the mapping solution. Sub- 
routine CLOSUR determines the integrated closure error of the solution. 


£a.ge_l 


ITERATION SUMMARY 


Heading 


Variable 


Description 


ITERATION 

SCALED MAXIMUM ERROR 

None 

None 

None 


V 

^m 

^CL 

ll 

Zo 


Iteration count 

Maximum scaled ^ror for the iteration 
Closure error \ ~ ^ 2 ^ 

Integration path #1 endpoint 
Integration path #2 endpoint 
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I*Iapped Duct Coordinate Page (2) 


Printed by MDAVIS 

Calculated by MDAVIS 

Options None 

Description 

Once the mapping iteration has terminated, the final solution and errors are 
printed for each wa.'l coordinate. 


Pa^e^ 1-2 I^P^ED DUCT_C00^2NATES 


Heading 


Variable 


Description 


TX 

TY 

X 

Y 

XC 

YC 

EX 

EY 

S 


X 

Y 

Xr 


"X 


Real part of t-plane pole location 
Imaginary part of t-plane pole location 
Image of T under mappping 

Input coordinates 

Error |X - X^. | 

Error [Y - Y^ | 

Arc length 
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Mesh Generation Page (1) 


Printed by COORD 

Calculated by COORD 

Options None 

Description 

The t-plane uniform mesh that is used 


for the coordinate generation is described. 


^ajge_l 


MESH PARAMETERS 


Heading 


Variable 


Description 


None 

None 

None 


JL 

DSTEP 

KN 


Number of streamwise steps 

Streamwise step size 

Number of uniform normal steps 
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10.4 Diagnostics and Failure Modes 


Numerous checks are performed during the course of the calculation. If a 
non-fatal or correctable error occurs a DIAGNOSTIC message is printed and the cal- 
culation continues. If a fatal error occurs a FAILURE Mode error is printed and 
the calculation stops. A DIAGNOSTIC message is printed of the form: 

** DIANGOSTIC NO. XX for 2-D COORDINATE OPERATOR and a FAILURE Mode 
message is of the form: 

** FAILURE NO. XX for 2-D COORDINATE GENERATOR where XX! refers to one 
of the conditions listed below. 

DIAGNOSTICS 

1) NUMERICAL SOLUTION OF SCHWARTZ-CRISTOFFEL TRANSFORMATION FAILED TO 
CONVERGE 

This error is detected in subroutine MDAVIS. It indicates that the scaled 
maximum error in the computed wall coordinates is greater than the input 
value ECONV after MAXIT iterations have been completed. By examining the 
ITERATION SUMMARY printed above the diagnostic message, one of three courses 
of action may be determined. 

a) The Scaled Maximum Error (SME) appears to be converging. Reset MAXIT 
and rerun the case. 

b) The SME has converged to a value different than zero. This can often 
be remedied by increasing the number of sub-steps (NSD) , employed in 
the normal direction integration. If this does not solve the problem, 
more wall definition coordinates may be needed. 

c) The SME is not converging. This often indicates that a poor initial 
potential flow solution was generated. 

2) UPPER AND LOWER WALLS NOT PARALLEL AT INLET. UPPER WALL FORCED PARALLEL, 

TO LOWER WALL 

This error is detected in subroutine ROTATM. It implies that the inlet upper 
and lower wall angles with respect to the horizontal differed by less than 
ten degrees but greater than 1.0-10. The upper wall endpoint is moved to 
force the walls parallel. 

FAILURE MODES 

1) MESH DISTORTION PARAMETER EQUALS = XXXXX 

This error is detected in subroutine ROBRTS. 
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2) MESH DISTORTION PARAMETER EQUALS = XXXXX 


This error is detected in subroutine DROBRT. 

3) LOWER WALL ANGLE = XXXXX DEGREES 
UPPER WALL ANGLE = XXXXX DEGREES 

WILL NOT FORCE PARALLEL IF DIFFERENCE IS > 10 DEG. 

This error is detected in subroutine ROTATM. It indicates that the inlet 
walls are not sufficiently parallel, and the program will not force the 
walls parallel to avoid drastically changing the geometry. 

4) DEGENERATE DERIVATIVE MAPPING FOR I = XXXXX 

dz 

This error is detected in subroutine MDAVIS. It implies that | < l.D-8 
at wall point XXXXX. 

5) INCONSISTENT OR INVALID INPUT 

This error is detected in main program CODUCT. 

Check input data set. 

6) READ ERROR ENCOUNTERED IN SUBROUTINE CORINP 

This error, detected while reading duct wall coordinates, indicates an error 
in the input data set. 

7) NUMBER OF INPUT POINTS EXCEEDS MAXIMUM (XXX) 

This error, detected in subroutine CORINP, indicates that too many wall values 
are defined. 

8) UNABLE TO COMPLETE APPROXIMATE SOLUTION 

This error is detected in subroutine ESTCOR, It implies that more than ten 
attempts have been made to compute a single potential line and is usually 
due to very large wall curvatures. 

9) INDEPENDENT STEP SIZE TOO SMALL 

This error is detected in subroutine DERIV3. It indicates that two con- 
secutively numbered wall coordinates are equal. 

10) IMSL LIBRARY ICSVKU FAILURE NO. XXX 

This error is detected in subroutine SMOOTH. It implies that the IMSL pro- 
gram ICSVKU cannot solve the spline problem. See IMSL manual to determine 
remedy. 
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11.0 GLOBAL STRUCTURE OF C0DUCT CODE 


This section of the manual is intended for the special user who wishes to 
modify the C0DUCT code or adapt to a different computer. The section provides a 
global overview of the code in terms of the principal tasks. These tasks are 
clearly labeled in the main program C0DUCT and agrees with the tasks listed on 
the Global Tree Structure Chart in Section 11.2. The global variables in labeled 
C0MM0N blocks are given in Section 11.3. Only the variables unique to the C0DUCT 
code are listed. Variables that are used by both the C0DUCT code and ADD code are 
listed in Section 6.0. Special problems associated with machine specific code are 
similar to those in the ADD code and are treated in Section 5.0. 
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11.1 Main Program CODUCT 


Object 

Main program for coordinate generator. 


Options 


IXFG 

= 

0 

Full solution 


= 

1 

Interpolate only 

lESTIM 

= 

0 

Approximate solution only 


= 

1 

Schwartz-Cristof f el transformation 

ICORD 

= 

0 

No coordinate output 


= 

1 

Coordinate output to disk files. 


Theory 

The control program CODUCT first calls subroutine RCNTRL (See Table 1) to read 
the user- input control parameters and options. Subroutine CORINP is then called to 
read the duct wall coordinates and, if requested, .will smooth the wall coordinate 
using a cubic-spline fitting algorithm. Subroutine ESTCOR is called to geometrically 
determine the approximate potential flow solution necessary to start the conformed 
mapping iteration procedure. Then Davis' algorithm to compute the Schwartz-Cristof f el 
transformation is invoked by calling subroutine MDAVIS. Finally, subroutine COORD 
is called to generate and output the coordinate mesh parameters to disk file(s). 
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11.2 Global Tree Structure By Task 



Read 

Control Input 


CODUCT 

RCNTRL 





Read 

and Smooth Duct Wall Data 


CODUCT 

CORINP 

SMARCL 

ARCLl 

SMOOTH 

ICSVKU 

Calculate Approximate Potential Flow Solution 

CODUCT 

ESTIMP 

ARCLl 

KURVTR 

DER1V3 




ESTCOR 

INSECT 

CROSSl 

CROSSl 



WALLV 

DER1V3 

SUNBAR 

UNBAR 



Calculate Sf'hwartz 

-Cristpffel Transformation 

CODUCT 

MDAVIS 

ROTATM 

INTNOR 

STEP 

STEP 




TTUP 

STEP 




CLOSUR 

INTNOR 

STEP 




INTSTR 

STEP 
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Calculate Coordinates and Metrics 

CQDUCT 

COORD 

CONSTR 

CDS 

ROBRTS 

DROBRT 






COORID 

NORLIN 

STEP 






STRSTP 

STEP 






CDVDN 

CORSTR 

Q2INTP 

BLKWRT 

NTRAN^ 






QPSTOR 

QPCURV 

INSECT 

CROSSl 



COORMD 

STRSTP 

STEP 






CDVDN 

CORSTR 

QPSTOR 

Q2INTP 

BLKWRT 

NTRAN^ 




1 


QPSTOR 

QPCURV 

INSECT 

CROSSl 

^ 

Transform from Uniform to Non-uniform niesh 

> 

CODUCT 

XFGRID 

INITQ 

CONSTR 

BLKRED 

NTRAN^ 





BLKRED 

NTRAN^ 






Q2INTP 

QPSTOR 

QPCURV 

INSECT 

CROSSl 




BLKWRT 

NTRAN^ 
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11.3 List of Labeled C0MM0N Blocks 


Name 


Object 


BSM0TH 

CESTP 

C00RC0 

C00RT 

IPDAVS 

NSDAVS 

0PDAVS 

TITLE 


Variables for spline smoothing 
Variables for approximate solution 
Control options and parameters 
Complex coordinates and derivative 
Complex solution variables 
Constants and parameters for mapping 
Intermediate Schwartz-Christof f el variables 
Run Title 
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List of Variables in C0MM0N/BSM0TH/ 
Variables for Spline Fitting 


Name 

Symbol 

Length 

Type 

Description 

A 


NXK 

Pv*4 

Integration constant 

B 


NXK 

R*4 

Integration constant 

CK 


ICK,3 

R*4 

Spline coefficients 

WK 


IWK 

R*4 

Work area 

XK 


NXK 

R*4 

Knot locations 

YPP 


NXNTD2 

R*4 

Second derivative of input curve 
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List of Variables in C0MM0N/CESTP/ 
Variables for Approximate Solution 


Name 

Symbol 

Length 

Type 

Description 

HT 


NXNTD2 

R*4 

Approximate duct height 

IL.IU 


NXNTD2 


Index lower /upper walls 

KKL.KKU 


NXNTD2 

R*4 

Curvature lower/upper walls 

KMEAN 

K 

NXNTD2 

R*4 

Curvature of mean line 

SL,SU 

^L’^U 

NXNTD2 

R*4 

Arc length lower/upper walls 

SLI.SUI 

^LI’^Ul 

NXNTD2 

R*4 

Arc length lower/upper walls 

SMID 

■5 

NXNTD2 

R*4 

Arc length mean line 

TH 

6 

NXNTD2 

R*4 

Angle of mean line with x axis 

VL,VU 


NXNTD2 

R*4 

Velocity lower /upper walls 

VLI.VUI 

\r"ui 

NXNTD2 

R*4 

Velocity lower/upper walls 

VMEAN 

V 

NXNTD2 

R*4 

Velocity on mean line 

XL,YL 


NXNTD2 

R*4 

Input coordinate lower wall 

XU, YU 


NXNTD2 

R*4 

Input coordinate upper wall 

XLI.YLI 


NXNTD2 

R*4 

Coordinates lower wall 

XUI,YUI 

^r^ui 

NXNTD2 

R*4 

Coordinates upper wall 


Note: subscript I denotes intersection of approximate potential line 

with wall 
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List of Variables in C0MM0N/C00RC0/ 
Control Options and Parameters 


Name 

Symbol 

Length 

Type 

Description 

DBS 

(An/An)j^ 


R*4 

Ratio uniform/nonuniform grid at 
wall 

IGRID 



1*4 

Grid option 

INUNIT 



1*4 

Output unit for nonuniform grid 

IRUNIT 



1*4 

Read unit 

ISM00T 



1*4 

Smoothing option 

lUUNIT 



1*4 

Output unit for uniform grid 

IWUNIT 



1*4 

Print unit 

IXFG 



1*4 

Transfer grid option 

JL 



1*4 

Number of output streamwise stations 

JLPTS 



1*4 

Number of input stations 

JXK 



1*4 

Number of knots for spline 

KL 



1*4 

Number of nonuniform streamlines 

KN 



1*4 

Number uniform streamlines 

L0P 



1*4 

Mesh distortion option 

RADR 

r 

r 


R*4 

Reference length 

TTL 



R*4 

Maximum t-plane coordinate lower 
wall 

TTU 



R*4 

Maximum t-plane coordinate upper 
wall 
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List of Variables in C0MM0N/C00RT/ 
Complex Coordinates and Derivatives 


Name Symbol 

T 

Length 

Type 

Description 

DZDTJ 1 

f— ) 

^*'(j+l 

j (d+1 
'dz \ 

KLL 

C*16 

Mapping derivative 

DZDTJP 1 

1 

KLL 

C*16 

Mapping derivative 

DZDTP2 

KLL 

C*16 

Mapping derivative 

DZDTl 1 

'dz\ 

,dt/ 

KLL 

C*16 

Mapping derivative 

ZJ 

Z^ 

KLL 

C*16 

Coordinate in duct plane 

ZJP 

J+1 

Z 

KLL 

C*16 

Coordinate in duct plane 

ZP2 

J+2 

Z 

KLL 

C*16 

Coordinate in duct plane 

Z1 

1 

Z 

KLL 

C*16 

Coordinate in duct plane 
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List of Variables in C0MM0N/IPDAVS/ 
Complex Solution Variables 


Name 

Symbol 

Length 

Type 

Description 

B 

b"-l 

NXNT 

R*8 

Location of poles 

BNEW 

b^ 

NXNT 

R*8 

Location of new poles 

TG 

tV-1 

NXNT 

C*16 

t-plane wall coordinates 

TT 


NXNT 

C*16 

New t-plane wall coordinates 

Z 

Z 

NXNT 

C*16 

Z-plane calculated wall coordinates 

ZC 

Z 

c 

NXNT 

C*16 

Z-plane input wall coordinates 
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List of Variables in C0MM(5N/NSDAVS/ 
Constants and Parameters for Mapping 


Name 

Symbol 

Length 

Type 

DM 

d 


R*8 


m 



EC0VV 

e 


R*8 


c 



IC 

0. + l.i 


1*4 

IC0RD 



1*4 

lESTIM 



1*4 

IPRINT 



1*4 

MAXIT 



1*4 

N 



1*4 

NBE 



1*4 

NCL,NCU 



1*4 

NLF 



1*4 

NMl 



1*4 

NNS 



1*4 

NUI 



1*4 

NIP0T 



1*4 

0NE 

1. + O.i 


C*16 

XM 

M 


C*16 

ZER 

0. + Oi 


C*16 


Description 

Automatic step size in approximate 
solution 

Convergence criteria 
Complex number 
Coordinate generator option 
Approximate potential flow option 
Print option 

Maximum number of iterations 
Number of wall points (NLF+2) 

Number of non-trivial angle changes 
Number of lower/upper wall elements 
Number of lower wall points 
N-1 

Number of additional steps in 
integration 

NLF+1 

Approximate potential flow option 
Complex 1.0 
Scale constant 
Complex 0.0 
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List of Variables in C0MM0N/TITLE/ 
Run Title 


Name 

S)nmbol 

Length 

Type 

Description 

ITITLE(I) 


18 

1*4 

Run title 
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List of Variables in C0MM0N/0PDAVS/ 


Name 

S3^bol 

Length 

Type 

Description 

ALPHA 

“j 

NXNT 

R*8 

Re(3j) + 1. 

BETAM 

Pj 

NXNT 

C*16 

Change in wall angle 

DELR 


NXNT 

R*8 

Difference in pole locations 

EXITA 

“e 

1 

R*8 

Duct exit divergence angle + n 

HEIGHT 

h 

1 

R*8 

Duct inlet height 

RATIO 


NXNT 

R*8 

Ratio of actual to calculated 
length 

THETAl 


1 

R*8 

Rotation of duct from real axis 

ZET 


NXNT 

C*16 

Estimated pole location 
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12.0 DETAILED DESCRIPTION OF C0DUCT CODE 


This section contains an alphabetic list of subroutines and a detailed 
description of each subroutine. Subroutines that are used by both the ADD code 
and the CODUCT code are listed and described in Section 7.0. The description of 
the subroutines which follow have the same format. This format consists of the 
object or purpose of the subroutine, any options used by the subroutine, and a 
list of variables not in a C0MM0N block which are used by the subroutine. Vari- 
ables in C0MM0N blocks are listed in Section 6.1 or Section 11.3. Following the 
list of variables is a brief description of the analysis performed by the sub- 
routine. 
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Name 

ARCLl 

BLKRED 

CDVDN 

C0BLK 

C0DUCT 

CL0SUR 

C0NSTR 

C00RD 

C00RMD 

C00R1D 

C0RINP 

C0RSTR 

CR0SS1 

DERIV3 

DR0BRT 

ESTC0R 

ESTIMP 

INSECT 

1NTN0R 

INTSTR 


12.0 DETAIL DESCRIPTION OF C0DUCT CODE 
12.1 List of Subroutines 

Obj ect 

Calculate arc length of input curve 
See Section 7.2 

Calculate streamline curvature 
Block data (IBM version) 

Main program (see Section 11.1) 

Calculate closure error in mapping 
Store fixed data in Ql, Q2 arrays 
Calculate coordinates and metrics 
Calculate coordinates J = 2, JL 
Calculate coordinates J = 1 
Read coordinate data 
Store coordinates in Ql array 
See Section 7.2 

Calculate 3 point central difference derivative 
See Section 7.2 

Determine locationof approximate potential line 
Calculate approximate pole locations 
Determine intersection of potential line and wall 
Calculate end point of potential line 
Calculate end point of streamline 



12.1 List of Subroutines (Cont'd) 


Name 

KURVTR 

MDAVIS 

N0RLIN 

QPKURV 

QPST0R 

Q2INTP 

RCNTRL 

R0TATM 

SMARCL 

SM00TH 

STEP 

STRSTP 

SUNBAR 

TTUP 

UNBAR 

WALLV 

XFGRID 


Obj ect 

Calculate curvature of input curve 

Solve Schwartz-Christof f el mapping 

Calculate single potential line 

Interpolates curvatures at output location 

Store Q parameters in Ql, Q2 arrays 

Interpolate from uniform to non-uniform mesh 

Reads user input control parameters 

Calculate duct rotation and scaling 

Cubic spline smoothing on arc length 

See Section 7.2 

Davis' integration formula 

Integrate each streamline one step 

Store X, y data into interpolation table T 

Update upper wall upstream point 

Lagrange table interpolation 

Approximate potential flow wall velocity 

Interpolate uniform to nonuniform grid 
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12.2 Description of Subroutines 


Subroutine ARCLl (X,Y,NPT,S) 


Object 

Calculate arc length of input curve 
Options 
None 


Symbols 

NPT 

S(I) 

X(I),Y(I) 




Number of input points 
Arc length 

Coordinates of input curve 


Theory 

The arc length of a curve is given by 




( 1 ) 
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Subroutine CDVDN (NPT,DZDT) 


Object 

Calculate streamline curvature 
Options 


None 



List of Symbols 


DZDT 

dZ/dt 

, Complex derivative of mapping 

NPT 


, Number of points in DZDT 

Q1(7,K) 

3v/3n 

, Streamline curvature 


Theory 

The magnitude of the potential flow velocity is 

V =1 ^ 

dZ 

Then the streamline curvature is given by 

K =-dV/dn 


( 1 ) 

( 2 ) 


The curvature is obtained by numerical differentiation using subroutine DER.IV3 for 
K = 2, KL -1. At the wall the streamline curvature is given by the wall curvature 
obtained from the input data. 
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Subroutine CL0SUR 


Obi ect 


Calculate closure error in mapping 
Options 
None 


List of Symbols 


NLF 

Z 

ZCL 

TT 



Theory 


Number of wall points 
Duct plane coordinates 
Closure error 
t-plane coordinates 


The solution is integrated from tj^ to + i by two paths to close the 

polygon. The closure error is defined by 


^CL“ '^1 Z (t|yj|_p + i)gj 


( 1 ) 
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Subroutine C0BLK 


BLOCK DATA 


Object 

Defines default values for program control 
List of Symbols 


NXNT 

Maximum 

NXNTD2 

Maximum 

maximum 

1ST 

Maximum 

NVK 

Maximum 

IPOINT 

Logical 

IWUNIT 

Logical 


number of wall definition 

number of wall definition 
number of potential lines 

number of streamlines 

number of knots to use in 

unit number to read from 

unit number to write to 


points total 

points for each wall, also 
permissible 

cubic-spline fit to wall data 
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Subroutine C0NSTR 


Object 

Store fixed data in Ql, Q2 arrays 
Options 

IGRID = 0 Uniform mesh 

= 1 Nonuniform mesh 

= 2 Both meshes 


Input Symbols 


DDS 

(An/An)j^ 

Mesh distortion parameter 

BETA 

An 

Uniform transverse step size 

JL 


Number of streamwise stations 

KL 


Number of nonuniform streamlines 

KN 


Number of uniform streamlines 

L0P 


Mesh distortion option 

SAVG 


Average length of duct 

Output 

Symbols 


Ql 

Q2 


Uniform mesh block data 
Nonuniform mesh block data 

Theory 




C0NSTR is a general setup program that stores information into the Ql and Q2 
data blocks. These stored variables are ones that do not change for J =1, JL. 

Q1(5,K) = 1, 

Q1(16,K) = 0. 

Q1(17,K) = 0. 

Q1(18,K) = n 

Q2(5,K) = 1. 

Q2(16,K) = dn/dn 
Q2(17,K) = n 
Q2(18,K) = n 
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Subroutine C00RD 


Object 

Calculate coordinates and metrics 

Options 

IGRID = 0 
= 1 
= 2 

Input Signals 

DDS An /An 

BETA An 

JL 

JLPTS 
KL 
KN 
L0P 
NNS 

Output Symbols 

Q1(I,K) 

Q2(I,K) 

Theory 

Once the mapping solution has converged, the location of the poles are known 
and the solution can be obtained for any interior point by integrating the Schwartz- 
Christoffel transformation. The streamwise integration step AS is defined in sub- 
routine C0NSTR by 

AS = (max - min (tjj.tL) ) / ( JL- 1 ) 

and the normal integration step is defined by 

An = 1/ (KN-I) 


Uniform grid output to lUUNIT 
Nonuniform grid output to INUNIT 
Output both grids 


Mesh distortion parameter 

Uniform step size 

Number of streamwise stations 

Number of points on wall 

Number of nonuniform streamlines 

Number of uniform streamlines 

Mesh distortion option 

Number of steps on potential line 


Coordinate data uniform grid 
Coordinate data nonuniform grid 


Then integration of dZ/dt with n constant produces a streamline and integration 
with S constant produces a potential line. 
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Subroutine C00RD (Cont'd) 


To start the coordinate calculation, C00R1D is called to Integrate the first 
potential line at the duct inlet. We note that dz/dt calculated by subroutine STEP 
is evaluated at the point 


^K+l/2 " ^^K+l 


(3) 


such that the metric is given by 




(4) 


Vk = 



(5) 


The remainder of the computation grid is constructed by integrating all the stream- 
lines in the streamwise direction one step using subroutine C00RMD. Again we note 
that the derivative is evaluated at the mid point so that 


/dz\^ j_ /dz\^~''^^ 

Idt/ ■ 2 lldt/ Idt/ 


( 6 ) 


V 




(7) 


The integration is continued to J = JL + 1/2. 
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Subroutine C00RMD (J) 


Obj ect 

Calculate coordinates J = 2, JL 
Options 


IGRID = 0 Uniform grid output to lUUNIT 

= 1 Nonuniform grid output to INUNIT 


® 2 Output 

Input Symbols 


DSTEP 

AS 

DZDTJ 

/dzf 

\dt/ 

DZDTJl 

\dt) 

J 


KL 


KN 



J 

ZJ 

Z 

ZJl 

Z'^-1 


both grids 

Streamwise step size 
Derivative at J 
Derivative at J-1 
Streamwise index 

Number of nonuhiform streamlines 
Number of uniform streamlines 
Coordinate at J 
Coordinate at J-1 


Output Symbols 

Q1 

Q2 

Theory 


Coordinate data uniform grid 
Coordinate data nonuniform grid 


The derivatives of the metrics (3V/3n, 9V/3S) are obtained by the three point 
difference formula. Thus we have 


and 




Idt/, 


( 1 ) 



J 

K 


( 2 ) 


III- 46 



Subroutine C00RMD (Cont'd) 


Since dZ/dt are known at the mid points we have 



_J 

2AS 



+ 1/2 


+ V 


J -*-3/2 



( 3 ) 


The derivative 3V/9n is obtained using subroutine CDVDN and the remaining variables 
are defined and calculated for the KN points on the uniform mesh by calling sub- 
routine C0RSTR. Subroutine Q2INTP interpolates from the KN uniform mesh points to 
the KL nonuniform mesh points. 
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Subroutine C00R1D 


Ojbect 

Calculate coordinates at J = 1 


Options 


IGRID = 0 Uniform grid output to lUUNIT 

= 1 Nonuniform grid output to INUNIT 

= 2 Output both grids 


Input Symbols 

DSTEP AS 

KL 

KN 

NNS 


Streamwise step size 
Number of nonuniform streamlines 
Number of uniform streamlines 
Number of steps in n integration 


Output Symbols 


DZDT2 

Derivative 

at J 

DZDT3 

Derivative 

at J 

Q1 

Coordinate 

data 

Q2 

Coordinate 

data 

Z2 

Coordinate 

at J 

Z3 

Coordinate 

at J 


uniform grid 
nonuniform grid 


Theory 


The first potential line is calculated by calling subroutine N0RLIN. Then 
derivatives of the metrics (3V/9n, 3V/3S) are obtained by 3 point difference formula. 
Thus we have 



( 1 ) 



( 2 ) 


The streamlines are integrated to J = 4 using subroutine STRSTP. Then we have 



+ 4 


V, 


3/2 


+ V 


5/2 

K 



(3) 
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Subroutine C00R1D (Cont'd) 


The derivative 3V/3n is calculated using subroutine CDVDN and the remaining variables 
are defined and calculated for the KN points on the uniform grid by calling sub- 
routine C0RSTR. Subroutine Q2INTP interpolates from the KN uniform grid to the KL 
nonuniform grid . 
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Subroutine C0RINP 


Object 

Read coordinate data 

Options 

ISM00T = 0 
= 1 

Input Symbols 

JLPTS 
JXK 
NLF 

Output Symbols 

XL,YL 
XU, YU 
ZC 

RADR 
Theory 

The subroutine reads the wall data in card image form. If ISM00T = 1, a cubic 
spline smoothing routine SMARCL will produce a set of JLPTS data points for each 
wall. This subroutine also prints the smoothed and unsmoothed wall data. 


Do not smooth wall data 
Smooth wall data 


Number of smoothed wall data points 
Number of knots in spline smoothing 
Number of input upper/lower wall data points 




Lower wall data points 
Upper wall data points 
Complex coordinates of wall data points 
Reference radius XU(1) 
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Subroutine C0RSTR (J,ZI,DZDT) 


Obi ect 

Store coordinates in Q1 array 
Options 

None 

Input Symbols 
J 

DZDT ^ ^ 

ZI 

Output Sybmols 

Q1 

Theory 


The following coordinate data are calculated at J. 


Q1(1,K) 

= 

J 

Im (Zj,) 

= R 

(1) 

Q1(2,K) 

= 

Re (Z^) 

= Z 

(2) 

Q1(3,K) 

= 

Re (dZ/dt)^ 

= 3R/3n 

(3) 

Q1(4,K) 

= 

Im (dZ/dt)^ 
K 

= 3R/3S 

(4) 

Q1(6,K) 

= 

l/ldZ/dtl^ 

= V 

(5) 

Q1(9,K) 

= 

r^j ^ 
Jo V 

= X 

(6) 

Ql(lO.K) 

= 

f\ ^ 
Jo V 

= Y 

(7) 

QKll.K) 

= 

Q1(10,K)/Q1(10,KN) 


(8) 

Q1(12,K) 

= 

Ztt ^ 

Jo y 

= A 

(10) 

Q1(13,K) 

= 

2tt R 


(11) 

Q1(14,K) 

= 

2tt 3R/3n 


(12) 

Q1(15,K) 

= 

2tt 9R/3S 


(13) 


Streamwise station number 
Derivative at J 
Coordinates at J 

Coordinate data block 
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Subroutine DERIV3 (X,Y,NX,NPT,DYDX,D2YDX2) 


Object 

Calculate 3 point central difference derivative 
Options 
None 

Input Symbols 

NPT 
NX 
X,Y 

Output Symbols 

DYDX d^/dx First derivative 

D2YDX2 d y/dx^ Second derivative 

Theory 

The finite difference formula are given by: 

/dY 

\dx) - 

- Y‘^'-lltrlY^rY-' 

\axr (x'+'-x')'^ + r(x'-x’-')* 


Point at which to evaluate derivative 

Number of data points for X and Y 

Table of NX independent/dependent variables 


r = 



(3) 


If I = 1 or NPT a diagnostic is printed 
" INPUT P0INT XX OUT OF RANGE" 
and both derivatives are set to 1.0. 
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Subroutine DERIV3 (X,Y,NX,NPT,DYDX,D2YDX2) (Cont'd) 

If |X^‘*'^- X^l or [X^- X^”^l < lO"^^ a diagnostic is printed 
"INDEPENDENT VARIABLE STEP SIZE LT l.E-15" 
and both derivatives are set to 1.0. 
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Subroutine ESTC0R (SAVG, NPOT) 


Object 

Determine location of approximate potential line. 


Option 


NIP0T 

0 

Program determines number ( 

□f lines 

> 

1 

NIP0T lines are calculated 


Input Symbols 




NP0T 

N 

P 

Number of potential lines 


SAVG 

S 

Average duct length 


XL,YL 

Vl 

Lower wall coordinates 


XU, YU 

Vv 

Upper wall coordinates 


Output Symbols 




HT 

h 

Height of duct 


IU,IL 




TH 

e 

Angle of mean line 


XLI.YLI 

%i»\i 

Coordinates potential line 

lower wall 

XUI,YUI 


Coordinates potential line 

upper wall 


Theory 


The object is to determine NP0T approximate potential lines in the duct where 
NP0T = NLF/3 initially. The first potential line intersects the duct at 
where the complex notation is used. 


Z =X + iY 


( 1 ) 
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Subroutine ESTC0R (Cont'd) 


We then construct a mean line Zm, (See Fig. 1), which satisfies the following con- 
conditions: 


^m,j = AS 


( 2 ) 


U1,J 




(3) 


(^U1,J ~ ^Ll,j) * (^m,j 


(4) 


It was found that the set of equations, Eq. (2) through (4) do not have a unique 
solution. Therefore Eq, (3) was replaced by a minimum condition on D where: 


D = I- 




(5) 


The algorithm consists of finding an angle 6^ which minimizes D. Thus we have from 
Eq. (1), 


^m,J " ^ m, j-i 5 



+ i sin 9 ^ ] 


( 6 ) 


A straight line normal to the mean line, from Eq. (4), is defined by the point 
Z’^m^j and the point. 


^ + AS«jcos {Qj +7T/2] + i sin (9^ + TT /2)^ 


(7) 


The intersections of the line (Z,Zj^ j) with the duct wall j, Z^^. j) is deter- 

mined using subroutine INSECT. Then D is calculated and checked for a minimum. An 
iteration procedure determines the 6^ which minimizes D. 
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Subroutine ESTC0R (Cont'd) 


When the iteration has converged, a check is made to determine if the Z„ ^ 

m y J 

potential line crosses the potential line inside the duct. If it does, 

the distance along the mean line AS is increased 

AS = AS • 1.2 (8) 


and the algorithm is repeated starting with Eq. (6). A maximum step AS is fixed by 
some fraction of the duct height d^j^. Thus 

, as] (9) 


AS = min [d 


m 


•UI.J "^LI,J 
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Subroutine ESTIMP 


Obj ect 

Calculate approximate pole locations 
Options 
None 

Input Symbols 
N 

NLF 
ZC 

Output Symbols 

B b^ Pole locations in plane 

TT t^ Pole locations in t plane 

Theory 

The arc lengths Sy,Sj^ to each pole in the duct (Z) plane is determined using 
subroutine ARCLl. Subroutine ESTC0R calculates the location of the approximate po- 
tential line and subroutine WALLV calculates the approximate potential flow veloc- 
ities (metrics) at the pole locations in the Z plane. Then the pole locations are 
given by: 


^U1 ~ J '^UI 

(1) 

0 


/•^L I 


tti =4 Vli ds 

(2) 

bui= exp (-TTtyj) 

(3) 

b^i = -exp (- TT tLi) 

(4) 


Total number of poles (2*NLF) 
Number of poles on each wall 
Z Duct coordinates 


111-58 



Subroutine INSECT (Arg. List) 


Object 

Determine intermine intersection of potential line and wall. 
Options 


lEXTRP = 1 

+ 1 

Argument List 

(X1,Y1),(X2,Y2) 
X,Y Z 

NPT 

(XI, YI) 

10 

lERR 


Extend last line segment for intersection 
Do not extend last line segment. 

Points defining potential line 

Points defining wall curve 

Number of (X,Y) points 

Intersection point 

Lower index of intersection point 

Error flag = 0 intersection found 

= -1 no intersection found 


Theory 

The input coordinates (X,Y) are searched for an intersection with (XI, Yl), 
(X2,Y2) using subroutine CROSSl which determines if an intersection occurs between 

Z (I 0) < Zj < Z(I 0 + I) (1) 
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Subroutine INTN0R (Arg. List) 


Object 

Calculate end point of potential line 
Options 


None 




Argument List 




ZETO 

^0 

Starting location 

in ? plane 

ZO 

Z 

o 

Starting location 

in Z plane 

ZETU 

'u 

Final location in 

C plane 

ZU 


Final location in 

Z plane 

NNS 


Number of steps 



Theory 

The starting location in the t plane is given by 


and the step size is given by 


At = i • I./(NNS-|) 


Then 


NNS-I ( 

= if 

■ rto+Ati 


dz 


dt 


J=l C'to + At(j-l) ) 

where the term in the bracket is evaluated using subroutine STEP. 


( 1 ) 


( 2 ) 


(3) 
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Subroutine Il^TSTR (Arg. List) 


Object 

Calculate end point of streamline 
Options 

None 

Argument List 


ZETO 


Starting location in . Q plane 

ZO 

"o 

Starting location in Z plane 

ZETU 

% 

Final location in C plane 

ZU 


Final location in Z plane 

List of Symbols 

NLF 


Number of lower wall points 

TT 

4 

Pole locations in t plane 

Theory 

The starting 

point 

in the t plane is given by 




Define the streamline by 


tp= Im (to) 


Then 


NLF ( .U~l + ♦n • ' H7 / 
^U" ^ 1 /f ^ ^ I 


where the bracket is evaluated using subroutine STEP. 


( 1 ) 


( 2 ) 


( 3 ) 
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Subroutine KURVTR (Arg. List) 


Object 

Calculate curvature of input curve 
Options 
None 

Argument List 


X,Y 

S 

NPT 

KURV 


X,Y Coordinates of input curve 

S Arc length of input curve 

Number of points on curve 
< Curvature of input curve 


Theory 

The principal curvature of a curve is given by 


- ^ 
" ds ds dS ds2 


Eq. (1) is evaluated by 3 point finite difference formula 


( 1 ) 
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Subroutine MDAVIS 


Object 

Solve Schwartz-Christoff el Mapping 
Options 

None 

Input Symbols 


EC0N 

e 

c 

Convergence criteria 

NLF 


Number of points on wall 

TT 

t 

Initial pole location t plane 

ZC 

Z 

Duct wall coordinates 


c 


Output Symbols 

TT t Final pole location in t plane 

Z Z Final pole location in Z plane 

Theory 

The flow chart for this subroutine is shown on Fig. 2 and takes place in the 
following steps: 

^t^p_l_Ini_tial^z^tion 

a) Calculate rotation constant M 

b) Calculate duct exit divergence angle 

c) Calculate Schwartz-Christof fel pole angle 
^t^p_2_Integr^a_te_Trans^ormat^on 

a) Integrate Schwartz-Christof fel transformation along each wall 
with a guess for the b^’s in the c plane using subroutine STEP. 

b) Integrate Schwartz-Christoff el transformation along far upstream 
potential line with a guess for b^’s in c plane using subroutine 
N0RLIN. 
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Subroutine MDAVIS (Cont'd) 


^t^p_3_U£da^t le ^o_les^ 

a) Update poles on lower wall in t plane by ratio 
of arc lengths. 

b) Update first pole on upper wall using Step 2b. 

c) Update poles on upper wall in t plane by 
ratio of arc lengths. 

d) Calculate pole location bj|^'s in ^ plane. 
^t£p_4_Che£k_C£nve£g£n£e_ 

a) Calculate absolute error for all poles 

“ Ih - ^cil 

b) Check convergence 

max (e^) < e 
i c 

c) Calculate closure error using subroutine CL0SUR 

d) If not converged repeat Steps 2, 3, 4 
If converged return 
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Fig. 2. Flow Chart for Subsonic MDAVIS 


• 1-10-S3-1* 
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Subr out ine H0HLIU ( J j Z ) 


Obj ect 

Calculate single potential line 
Options 
None 

Input Symbols 


J 


Calculate Jth potential line 

KN 


Number of output stations 

NNS 


Number of integration steps 

TT 

t 

i 

Pole locations in t plane 

ZO 

Z 

0 

Initial Z location 

Output Symbols 



Z1 

Z 

K 

Coordinates of potential line 

DZDTl 

(dZ/dt) 

K 

Derivative 


Theory 

After a converged solution for the pole locations is obtained, this subroutine 
integrates the potential line at the Jth station in NNS steps and outputs the 
coordinates and derivatives of KN stations. Then let us choose 


NNS = NSD* (KN - I ) + K N 


■ ( 1 ) 


where NSD is the number of integration steps per output station. The integration 
starts at in the Z plane and t^ given by 

tg = DSTEP (j-|)-t|_/2+0+i (2) 


in the t plane. The parameter tj^ is chosen in the approximate coordinate calcula- 
tion to center the pole distributions about plus and minus values. 
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Subroutine N0RLIM (Cont'd) 


The integration step is then given by 


At =1/ (NNS-I)* i 


( 3 ) 


Then we have the recursion formula 


z. 



( 4 ) 


Z», = Z 


K-l 


I 


t^-,+At*L 


+At- {l-D 


dz 

dt 


dt 


( 5 ) 



( 6 ) 
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Subroutine QPKURV 


Ob.i ect 


Interpolates wall curvature at output location 


Options 

None 

Input Symbols 


KKL, KKU 


XL,YL 


XU, YU 


SL,SU 


Q1 


Output Symbols 


RHSK3) 

K^CJ) 

RTS1(3) 

Ky(j) 


Theory 


Curvature of lower/upper wall 
Input coordinates lower wall 
Input coordinates upper wall 

Arc length lower /upper wall 
Coordinate data 

Curvature lower wall 
Cuirvature upper wall 


The streamline curvature KKL and KKU is known at the input data points (XL,YL) 
and (XU, YU) respectively^ The coordinates are known at station J for equal stream- 
wise steps DSTEP. Let (Xj^.Y^) and (Xy,Y^) be the lower and upper wall coordinates 
at station J obtained from the Q1 array. A straight line is passed through these 
points and a serarch of the input coordinates is made using subroutine INSECT to 
determine the intersection on each wall. Subroutine INSECT returns an interpolation 
parameter which is used to calculate K^(J),Ky(J) 
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Subroutine QPST0R(J) 


Object 

Store Q parameters in Ql, Q2 arrays 
Options 


IGRID = 0 

= 1 

= 2 

IXFG = 0 

= 1 

Input Symbols 

DSTEP AS 

JL 

KL 

RADR r 

r 

Output S 3 nnbols 

RHS1,RMS1,RTS1 
RHS2,RMS2, RTS2 
QPARMl 
QPARM2 


Uniform grid 
Nonuniform grid 
Both grids 
Uniform grid 
Interpolate only 

Streamwise step size 
Number of potential lines 
Number of streamlines 
Reference radius 

Wall coordinate data uniform grid 
Wall coordinate data nonuniform grid 
Parameters for uniform grid 
Parameters for nonuniform grid 


Theory 

The wall, coordinate data and grla parameters are calculated by this subroutine. 
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Subroutine Q2INTP 


Object 


Interpolate from uniform grid to nonuniform grid 
Options 
None 

Input Sybmols 


KL 

Number of output 

streamlines 

KN 

Number of input 

streamlines 

Q1 

Coordinate data 

for uniform grid 


Output Symbols 

Q2 Coordinate data for nonuniform grid 

Theory 

The normal coordinate for a uniform grid is Q1(19,K) K=1,KN and the normal 
coordinate for the nonuniform grid is Q2(18,K) K=1,KL. The Q2 variables are ob- 
tained from the Q1 variables by linear interpolation using the normal coordinate 
as the independent variable. 
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Subroutine RCNTRL 


Object 

Reads user input control parameters 
Options 

None 

Input Signals 

See input data Section 10.2 
Theory 

This subroutine reads the input control parameters, checks for inconsistencies 
and prints the input data. 
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Subroutine R0TATM 


Object 

Calculate duct rotation and scaling 
Options 

None 

Input Symbols 

ZC Input duct coordinates 

Output Symbols 

M Rotational constant 

Theory 

Far upstream of the duct inlet ^ » and the Schwartz-Christof f el 

tlon reduces to 

dz _ M 
dC " C 

Integrating Eq, (1) we have 

Z = Min C + Zq 

The tranf ormation to the t plane is given by 

I n C 7T ( i - 1 ) 

and :Eq. (2) becomes a duct with' parallel walls. 

Z = M 7T (i-t) + Zq 

Subtracting the lower wall from the upper wall we have 

Z u ~ Z^= - M 7T i 


transf orma- 

( 1 ) 


( 2 ) 


(3) 


(4) 


(5) 
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Subroutine ROTATM (Cont'd) 


The height of the duct is given by 


H = 



( 6 ) 


Hence 




(7) 


where 6 is the angle of the duct with respect to the real axis. 
Then solving for M using Eq. (5) and Eq. (7) we have 



id 


( 8 ) 


Since this solution requires parallel walls, this subroutine will modify the Nth 
point of the data set to insure parallel walls at the inlet so that the correct duct 
height H can be determined at the inlet. 
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Subroutine SMARCL (Arg. List) 


Object 


Cubic spline smoothing 

on arc length 

Options 



None 



Argument List 



JX 


Number of input coordinate points 

JXB 


Number of ouptut coordinate points 

JXK 


Number of knots in spline fit 

X,Y 

X,Y 

Input coordinates 

S 

X 

Arc length along input curve 

SB 

S 

Arc length along output curve 

XB,YB 

X,Y 

Output smoothed coordinates 

Theory 




The arc length along the input curve is calculated using subroutine ARCLl. 
Then an increment of arc length is defined by 

AS = (S (JX) - S(l) ) / (JXB-1) (1) 


The £uryes X(S) and Y(S) are smo£thed using subroutine SM00TH which returns X(S) 
and Y(S) for JXB points spaced AS in length 
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Subroutine STEP(ZETD1,DZETD,DZ) 

Ob j ect 

Calculate Integration Step For Schwartz-Christof f el Transformation 
Options 


None 




Variables 




B(K) 

= 


Location of pole in ? plane 

BETM(K) 

= 


Turning angle in Z plane 

DZ 

= 

m 

Step size in Z plane 

DZETD 

= 


Step size in ^ plane 



m 

Step size in t plane 

XM 

= 

M , 

Scale factor 

ZETDl 

= 

m 

Initial ^ 

ZETD2 

= 

^m+l » 

Final r. 

NBE 

= 

N , 

Number of poles 

GAMA 

= 

''^m ’ 

Exit divergence angle 


Theory 

The second order integration formula evaluated at the mid point is given by 
Davis Ref. 1 as 

/ -Qj/TT+I -Oj/tT+I 

/dZ\ _ M .QM/tt n ) (^m + I ~^i) 

'^Wm+I/2 ^M+I/2 1=1 I + 1 ) 



where 




^ M +1 

= — (r + 

M + I/2 2 '^M + t 



(2) 

( 3 ) 
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Subroutine STEP (Cont'd) 


The transformation to the t plane is given by 




1/2 


Then we have 




where At is chosen by the input c's. 
in 




References 


1. Davis, R. T.: Numerical Methods for Coordinate Generation Based on 

Christoff el Transformations. 


( 4 ) 


(5) 


( 6 ) 


Schwartz- 
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Subroutine STRSTP 


Object 

Integrate each streamline one step 
Options 

None 

Argument List 


J 

ZI 

ZO 

DZDTO 


,J+1 


(dZ/dt) 


J+1/2 


Streamwise station 
Coordinates at J 
Coordinates at J+1 
Derivative at mid point 


Theory 

This subroutine integrates K = 1, KN streamlines one step 


Z 



■+ 


+ 1/2 

dt^K 


At 


using subroutine STEP. 



Subroutine SUNBAR (X,Y,T,NPT,NORDER) 


Object 

Stores X,Y data into interpolating Table T 
Options 


None 


Input Symbols 

X,Y x,y Independent variable arrays 

NPT Number of data pairs (X,Y) 

NORDER Interpolation order (=1, 2 or 3) 

Output Symbols 


T Output interpolation table . 

Theory 

The data is stored into T as follows 

T(l) = 1. 

T(2) = NORDER 
T(3) = NPT 
T(4) = 0. 

T(J + 4) = X(J), J = 1,NPT 
T(J + NPT + 4) = Y(J), J = l.NPT 
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Subroutine TTUP (ITER,ZU) 


Object 

Update upper wall upstream point 
Options 


None 




Input Symbols 




ITER 

V 

Iteration number 


THETAl 

®1 

Angle of duct rotation 

from real axis 

TT 

t? 

X 

Location of poles in t 

plane 

ZC 

Z . 

Cl 

Input wall coordinates 


zu 

"u 

End point of potential 

line integration 

Output Symbols 




TT(N) 

v+1 

Updated t plane coordinate at point N 


Theory 

This subroutine updates the corner point Z to close the polygon by jumping 
from the lower wall to the upper wall as shown’^^in Fig. 3. With known t., the point 
Z^ is determined by integrating along the path (Z^ 2 _»A,Z^)» The error in closing the 
polygon is given by 





( 1 ) 


\H-1 

The update tjl^ 
points tj^ and t 


N 


is determined in the following manner, 
given by 


t 


I 

I 



Let us define upstream 


( 2 ) 


t N = /| + i 


(3) 
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Subroutine TTUP (Cont'd) 


where a is a parameter chosen to move t^^ sufficiently upstream to approximate the 
limiting case t -> - <» 


Z = 7tM (i-t) + Zp 


(4) 


Then 


z'^-z' =-7TMi 


(5) 


* T t 

The point is determined by integrating along the path Z^) and is deter- 

rained from Eq. (5). A ratio of wall lengths is defined 



and the update for t^^ is given by 


I 


-•k+RW-C) 


( 6 ) 


(7) 


\ 
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2 PLANE (Z = X + iY) 



X 


t PLANE (t = s^-in) 



Fig, 3. Update for Corner Point 


81 - 10 - 63-3 
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Subroutine UNBAR (T,1K,XIN,Y1N,ZZ,KK) 


Object 

Interpolate a univariate or bivariate table. 


Input Symbols 
T 


Name of the array which contains the table values. 


IK 


Element of the array at which the table starts. If 
you have only one table in the array, IK=ONE. 


XIN 


Independent variable in the X-sense. 


YIN 


Independent variable in the Y-sense. If the table 
is a univariate, then YIN is zero. 


Output Symbols 

ZZ 

KK 


Dependent variable 
Off Table indicator 

0 Normal evaluation 

1 Off On X Min. 

2 Off On X Max. 

3 Off On Y Min. 

4 Off On X Min. and Y Min. 

5 Off On X Max. and Y Min. 

6 Off On Y Max. 

7 Off On X Min. and Y Max. 

8 Off On X Max. and Y Max. 

Less Than 0, Table set up wrong. 


Theory 

If either variable is off the table, UNBAR will return the corner value.. 
This implies that UNBAR will not extrapolate and does not recognize any discon- 
tinuities. The table must be set up as follows-all numbers are in floating point 
mode. 


T(IK) 
T(IK+1) = 
T(IK+2) = 
T(IK+3> = 
T(IK-H-) = 
T(IK++) = 
T(IK++) = 


Curve No. 

Degree of Interpolating (1, 2, 3) 

NX. No. of X values 

NY. No. of Y values, (in univariate make zero) 

X values in ascending order. 

Y values in ascending order. 

Z values. Put them in following order- (Z(l.l) ,Z(1,2) , 

Z(l,3) Z(1,NY),Z(2,1),Z(2,2) Z(2,NY) Z(NX,1), 

Z(NX,2) Z(NX,NY). For bivariate only. 
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Subroutine UNBAR (Cont'd) 


A Lagrongian interpolation polynomial of degree 1, 2 or 3 will be uSed for the 
interpolation depending upon T(IK+1). 
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Subroutine WALLV 


Object 

Approximate potential flow velocity 
Options 


None 

Input Symbols 

HT 

h 

9 

Approximate duct height 

IL,IU 


» 

Indices for lower/upper wall potential 
line 

NP0T 

N 

P 

f 

Number of potential lines 

SAVG 

? 

9 

Average duct length 

SL,SU 


9 

Arc length lower /upper wall 

TH 

0 

9 

Mean line angle 

XL,YL 


9 

Coordinates lower wall 

XU, YU 

^U’^u 

9 

Coordinates upper wall 

XLI.YLI 


9 

ESTC0R coordinates lower wall 

XUI.YUI 

^ir^ui 

9 

ESTC0R coordinates upper wall 

Output Symbols 

VL.VU 

Vt ,V 
L’ U 

9 

Velocity lower/upper wall 


Theory 

For each computed potential line, the wall velocity can be estimated as 
follows 


V=l/n (1) 


and the curvature of the mean line by 


K = 


dS 


( 2 ) 
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Subroutine WALLV (Cont'd) 


Then define 


l-K/(2V) 
l + K/(2V) 


and the approximate velocity at the wall is given by 


( 3 ) 


V 


UI 


2 

I + <f> 


V 


(4) 




(5) 


Arc lengths Slx and may be calculated using subroutine ARCL for the lower and 

uppwer walls defined by and linear interpolation is used 

with arc length as an independent variable to interpolate V and V„ from the table 
„ Lb 
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Subroutine XFGRID 


Object 


Interpolate uniform to nonuniform grid 
Options 


None 

Input Symbols 

DBS 

KL 

Q1 

Output Symbols 

Q2 


Mesh distortion parameter 
Number of output streamlines 
Input coordinate data array 


Output coordinate data array 


Theory 

If IXFG option is turned on, only this subroutine is called by the main program 
C0DUCT. This subroutine reads in input coordinate file from unit lUUNIT, interpolates 
the Q1 array with KN uniformly spaced streamlines to obtain the Q2 array with KL 
nonuniformly spaced streamlines and stores the output on unit INUNIT. 
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